Análisis de la pelea final de Los Vengadores: Endgame.

Introducción

Para la realización de este ejercicio he empleado la base de datos avengers.csv de la librería flexplot. He modificado la base de datos para solo tener en cuenta las variables que se emplean en este video y les he cambiado el nombre. Si queréis la base de datos original, la podeís encontrar en este enlace.

Esta base de datos contiene originalmente los atributos de combate de 812 luchadores en la pelea final de Los Vengadores: Endgame. Para este ejercicio solo se emplearán 3 atributos, que son los siguientes: minutes.fighting (renombrado a T_Batalla), son los minutos que un luchador es capaz de durar en la batalla hasta que muere o se rinde; willpower (renombrado a Determinacion), que según la documentación, está determinado por el periódo de tiempo que el luchador es capaz de esperar en la DMV (aquí en España es la DGT) para que le entreguen el carnet de conducir; y finalmente, injuries (renombrado a Heridas) que es el número de heridas que pueden soportar.

Escena de la batalla final de Los Vengadores: Endgame

Escena de la batalla final de Los Vengadores: Endgame

En este ejercicio vamos a tratar de modelizar el tiempo en batalla (T_batalla) en función de las variables Determinación y Heridas empleando la regresión gamma.

En primer lugar vamos a cargar el banco de datos, y siguiendo la forma en que trata a la variable Heridas en el vídeo del enlace del primer párrafo, vamos a categorizar en tres grupos la variable Heridas mediante la función cut().

datos <- read.csv("avengers_mod.csv");datos <- datos[-437,] ## valor negativo de determinación
datos$Heridas <- cut(datos$Heridas, breaks = c(-1,2,4,5),labels = c("0-2","3-4","5+"))

Con esto realizado, podemos empezar con la parte de estadística descriptiva del banco de datos.

Análisis previo descriptivo del banco de datos

En primer lugar vamos a obtener un resumen numérico de las variables mediante la función summary().

summary(datos)
##    T_batalla     Determinacion    Heridas  
##  Min.   : 1.10   Min.   :  4.00   0-2:206  
##  1st Qu.: 8.70   1st Qu.: 48.00   3-4:374  
##  Median :10.40   Median : 60.00   5+ :231  
##  Mean   :12.38   Mean   : 60.07            
##  3rd Qu.:13.15   3rd Qu.: 72.00            
##  Max.   :93.70   Max.   :115.00

Observamos que no hay ningún valor ausente en ninguna de las variables. La mayoría de los valores de nuestra variable respuesta T_batalla se concentran en el intervalo [1.10,13.12] quedando algunos de estos valores muy alejados del resto. En cuanto a los valores de la variable Determianción podemos intuir que son bastante simétricos, con media en 59.99. Finalmente, en cuanto a la variable Heridas podemos decir que la hay un número similar de luchadores con 0-2 heridas, que con 3-4 y que con 5+,siendo el grupo con 3-4 heridas el más numeroso.

Pasamos ahora a la parte gráfica del análisis descriptivo. Vamos a comenzar esta parte obteniendo el histograma T_batalla.

ggplot(data = datos, aes(T_batalla)) +
  geom_histogram(color = "dodgerblue", fill = "skyblue") + 
  labs(title = "Figura 1. Histograma de T_batalla") +
  theme_bw()

El histograma cumple con una de las particularidades de la distribución gamma: es asimétrica por la derecha. Observamos que las probabilidades de que un luchador esté mucho tiempo peleando son muy bajas, durando la mayoría de ellos tiempos comprendidos entre los 0 y los 25 minutos.

Ahora vamos a ver qué tipo de relación hay entre las variables. Para ello, vamos a obtener el diagrama de dispersión entre T_batalla y Determinacion .

ggplot(data = datos, aes(x = Determinacion, y = T_batalla)) +
  labs(title = "Figura 2. Diagrama de dispersión T_batalla frente a Determinacion") +
  geom_point() + theme_bw()

A través de este diagrama podemos intuir que cuanto mayor es la Determinacion mayor es el tiempo en batalla. Sin embargo, no está del todo claro. Así que vamos a ver si hay alguna relación si añadimos la variable Heridas y obtenemos un diagrama de dispersión para cada grupo.

ggplot(data = datos, aes(x = Determinacion, y = T_batalla)) +
  labs(title = "Figura 3. Diagramas de dispersión de T_batalla frente a Determinacion por grupos de Heridas") +
  geom_point() + facet_wrap(~Heridas)  + theme_bw() +
  theme(plot.title = element_text(size=11))

Ahora la relación parece bastante más clara. A mayor Determinación máyor es T_batalla, solo que la variable Heridas disminuye el efecto que tiene la Determinacion, de modo que a mayor número de heridas, menos importante es el efecto de Determinacion. Es decir, parece que hay interacción entre las variables Determinacion y Heridas.

Modelización Lineal

Aunque sepamos que estamos analizando una variable continua no negativa, ajustamos los datos a un modelo lineal para observar si se ajusta adecuadamente.

mod.lineal <- lm(T_batalla ~ ., data=datos)
summary(mod.lineal)
## 
## Call:
## lm(formula = T_batalla ~ ., data = datos)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -10.808  -2.913  -0.711   1.286  75.261 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   11.49952    1.07994  10.648  < 2e-16 ***
## Determinacion  0.10676    0.01469   7.269 8.55e-13 ***
## Heridas3-4    -6.33453    0.64696  -9.791  < 2e-16 ***
## Heridas5+     -9.15367    0.71184 -12.859  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.376 on 807 degrees of freedom
## Multiple R-squared:  0.2411, Adjusted R-squared:  0.2383 
## F-statistic: 85.48 on 3 and 807 DF,  p-value: < 2.2e-16
par(mfrow=c(2,2)); plot(mod.lineal);par(mfrow=c(1,1))

Como podemos observar en los plots de diagnóstico, los residuos no se ajustan adecuadamente. En primer lugar, el gráfico residuals vs. fitted no sugiere una relación lineal perfecta, puesto que hay residuos que no siguen la línea horizontal. muchos de los residuos se distribuyen a lo largo de la línea horizontal. Por otra parte, el QQ plot muestra como la mayoría de los residuos caen en la línea de referencia, pero no así el extremo superior. Esto sugiere que los datos siguen una distribución de cola larga. Tampoco parece que se cumpla adecuadamente la hipótesis de homocedasticidad. El gráfico del leverage muestra también muchos valores extremos. Por tanto, la regresión gaussiana no parece la mejor alternativa para modelizar estos datos.

Modelización con Regresión Gamma

Como estamos analizando una variable continua no negativa en la que de forma general su varianza parece aumentar con la media, parece razonable emplear un ajuste con una regresión gamma.

La variable respuesta \(Y_{ij}\) representa el número de minutos que aguanta en combate el luchador i perteneciente al grupo j de la variable Heridas.

La variable respuesta \(Y_{ij}\) se distribuye Gamma de parámetros \(\nu\) y \(\lambda\):

\[ Y_{ij} \sim Ga(\nu_{ij}, \lambda_{ij}) \]

La componente sistemática del modelo está formada por 2 variables explicativas (Determinación, variable cuantitativa continua, y Heridas, variable categórica ordinal).

Sabemos que \(\lambda = \nu / \mu\), por lo que \(E(Y_{ij})=\frac{\nu}{\lambda}=\mu_{ij}\). Para relacionar el predictor lineal con la respuesta media, podemos emplear 3 funciones link:

  • La función inversa: \(g(\mu_{ij}) = \frac{1}{\mu_{ij}}\)
  • La función log: \(g(\mu_{ij}) = \log \mu_{ij}\)
  • La función identidad: \(g(\mu_{ij}) = \mu_{ij}\)

De tal forma que el predictor lineal queda: \[ g(\mu_{ij}) = \beta_0 + \beta_1 x_{i} + \gamma_1 d_j +\gamma_2 d_j \]

donde \(x_{ij}\) es la Determinación del individuo i y dj es una variable indicadora para la que \(\gamma_1 \cdot 1\) si Heridas = 3-4 y \(\gamma_2 \cdot1\) si Heridas = 5+.

Comparación de modelos: Validación del ajuste y Capacidad predictiva

A continuación se muestra una tabla comparativa a nivel ajuste y predicción de los diferente modelos.

Link Modelo %D explicada AIC RMSE
Inverso Completo 44.67 4549.9 53.6
Logaritmo Interacción 45.28 4544.8 53.6
Identidad Interacción 45.22 4545.7 53.5

Por lo tanto, podemos determinar que no hay ningún modelo mucho mejor que otro en términos de ajuste y predicción, siendo el modelo con la función link logarítmica la que presenta, aunque por poco, un valor de AIC y DEVIANCE menor. Además, es el que mejor interpretabilidad tiene.

Finalmente, podemos ver graficamente los valores ajustados y observados del modelo seleccionado.

par(mfrow=c(1,1))
plot(log.interac$fitted.values, datos$T_batalla, main='Mejor modelo', xlab='Valores Ajustados', ylab='Valores Observados')
abline(0,1)

Modelo final

Por ello, elegimos este modelo como modelo final. \[ Y_{ij} \sim Ga(\nu_{ij}, \lambda_{ij}) \] \[ log(\mu_{ij}) = 2.06 + (0.012 -0.007d_{(3-4)}-0.005d_{(5+)})x_{i} + 0.03d_{(3-4)} -0.33 d_{(5+)} \]

con \(X_i\) como variable Determinacion y \(d_j\) como variable indicadora de Heridas.

Así pues: \[ \widehat{E(Y_{ij})} = \mu_{ij} = exp(2.06 + (0.012 -0.007d_{(3-4)}-0.005d_{(5+)})x_{i} + 0.03d_{(3-4)} -0.33 d_{(5+)}) \] Por lo tanto, el modelo se interpreta como cambios porcentuales en \(E(Y_{ij})\) por cada incremento de una unidad en \(X_i\) (% = \(100 \cdot (e^{\beta_i})\)).

Análisis del BMI mediante regresión gamma

Introducción

El siguiente análisis consiste en la modelización mediante regresión gamma de los datos de la American Time Use Survey (ATUS) Eating & Health Module de 2014. Concretamente, se modelizará el BMI en función del resto de variables.

Se puede encontrar más información sobre los datos en su pagína web, aunque el dataset empleado en este caso es una conversión a CSV que se puede encontrar en Kaggle.

En primer lugar, cargamos los datos (eliminando aquellas observaciones para las que el BMI reportado sea 0, ya que es imposible):

bmi_data <- read_csv("data/ehresp_2014.csv", show_col_types = FALSE) %>%
filter(erbmi > 0)

Análisis descriptivo

bmi <- bmi_data
str(bmi, give.attr=FALSE)
## spec_tbl_df [10,637 × 37] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
##  $ tucaseid   : num [1:10637] 2.01e+13 2.01e+13 2.01e+13 2.01e+13 2.01e+13 ...
##  $ tulineno   : num [1:10637] 1 1 1 1 1 1 1 1 1 1 ...
##  $ eeincome1  : num [1:10637] -2 1 2 2 1 1 1 1 1 1 ...
##  $ erbmi      : num [1:10637] 33.2 22.7 49.4 31 30.7 ...
##  $ erhhch     : num [1:10637] 1 3 3 3 3 1 3 3 3 3 ...
##  $ erincome   : num [1:10637] -1 1 5 5 1 1 1 1 1 1 ...
##  $ erspemch   : num [1:10637] -1 -1 -1 -1 1 5 -1 -1 5 -1 ...
##  $ ertpreat   : num [1:10637] 30 45 60 65 20 30 30 117 80 35 ...
##  $ ertseat    : num [1:10637] 2 14 0 0 10 5 5 10 0 20 ...
##  $ ethgt      : num [1:10637] 0 0 0 0 0 0 0 0 0 0 ...
##  $ etwgt      : num [1:10637] 0 0 0 0 0 0 0 0 0 0 ...
##  $ eudietsoda : num [1:10637] -1 -1 -1 -1 1 -1 -1 -1 2 1 ...
##  $ eudrink    : num [1:10637] 2 2 1 1 1 1 2 2 1 1 ...
##  $ eueat      : num [1:10637] 1 1 2 2 1 1 1 1 2 1 ...
##  $ euexercise : num [1:10637] 2 2 2 1 1 2 1 1 2 2 ...
##  $ euexfreq   : num [1:10637] -1 -1 -1 5 2 -1 3 6 -1 -1 ...
##  $ eufastfd   : num [1:10637] 2 1 2 2 1 1 1 2 1 1 ...
##  $ eufastfdfrq: num [1:10637] -1 1 -1 -1 3 3 1 -1 2 5 ...
##  $ euffyday   : num [1:10637] -1 2 -1 -1 1 2 2 -1 1 2 ...
##  $ eufdsit    : num [1:10637] 1 1 1 1 1 1 1 1 1 1 ...
##  $ eufinlwgt  : num [1:10637] 5202086 29400000 26000000 17500000 3661280 ...
##  $ eusnap     : num [1:10637] 1 2 2 1 2 2 2 2 2 2 ...
##  $ eugenhth   : num [1:10637] 1 2 5 4 3 2 2 3 1 4 ...
##  $ eugroshp   : num [1:10637] 1 3 2 1 2 3 1 1 1 1 ...
##  $ euhgt      : num [1:10637] 60 63 62 69 71 65 63 70 65 70 ...
##  $ euinclvl   : num [1:10637] 5 5 5 5 5 5 5 5 5 5 ...
##  $ euincome2  : num [1:10637] -2 -1 2 2 -1 -1 -1 -1 -1 -1 ...
##  $ eumeat     : num [1:10637] 1 1 -1 1 -1 1 1 1 1 -1 ...
##  $ eumilk     : num [1:10637] 2 2 -1 2 -1 2 2 2 2 -1 ...
##  $ euprpmel   : num [1:10637] 1 1 2 1 2 3 1 1 1 2 ...
##  $ eusoda     : num [1:10637] -1 -1 2 2 1 2 -1 -1 1 1 ...
##  $ eustores   : num [1:10637] 2 1 -1 1 -1 2 1 1 3 1 ...
##  $ eustreason : num [1:10637] 1 2 -1 1 -1 5 3 4 1 1 ...
##  $ eutherm    : num [1:10637] 2 2 -1 2 -1 2 2 2 2 -1 ...
##  $ euwgt      : num [1:10637] 170 128 270 210 220 200 155 180 170 282 ...
##  $ euwic      : num [1:10637] 1 2 2 1 2 2 -1 -1 -1 -1 ...
##  $ exincome1  : num [1:10637] 2 0 12 0 0 0 0 0 0 0 ...

Tenemos un banco de datos con 10637 observaciones y 37 variables. Las dos primeras variables son variables de identificación, y la última es de etiquetado, por lo que las eliminamos por simplicidad.

bmi <- bmi[, -c(1, 2, 37)]

Respecto a las variables restantes, podemos destacar que:

  • “erbmi” es una variable cuantitativa continua, y constituye nuestra variable respuesta.
  • “ertpreat”, “ertseat”, “euexfreq”, “eufastfdfrq”, “eufinlwgt”, “euhgt”, “euwgt” son variables cuantitativas.
  • El resto de variables son cualitativas, por lo que han de ser convertidas a factor.

En primer lugar, observamos las correlaciones entre las variables.

corrplot(cor(bmi), type = "upper")

Se aprecian correlaciones entre las variables que aportan información similar. En el caso del BMI una correlación importante con el peso. Respecto a las variables cuantitativas:

no.factor <- c("erbmi", "ertpreat", "ertseat", "euexfreq", "eufastfdfrq", 
               "eufinlwgt", "euhgt", "euwgt")

summary(bmi[,no.factor])
##      erbmi          ertpreat         ertseat         euexfreq    
##  Min.   :13.00   Min.   :  0.00   Min.   : -2.0   Min.   :-2.00  
##  1st Qu.:23.60   1st Qu.: 30.00   1st Qu.:  0.0   1st Qu.:-1.00  
##  Median :26.60   Median : 60.00   Median :  4.0   Median : 2.00  
##  Mean   :27.77   Mean   : 65.85   Mean   : 16.9   Mean   : 2.27  
##  3rd Qu.:30.70   3rd Qu.: 90.00   3rd Qu.: 15.0   3rd Qu.: 4.00  
##  Max.   :73.60   Max.   :508.00   Max.   :990.0   Max.   :38.00  
##   eufastfdfrq       eufinlwgt             euhgt           euwgt      
##  Min.   :-2.000   Min.   :   756844   Min.   :56.00   Min.   : 98.0  
##  1st Qu.:-1.000   1st Qu.:  3511427   1st Qu.:64.00   1st Qu.:145.0  
##  Median : 1.000   Median :  6023405   Median :66.00   Median :170.0  
##  Mean   : 1.158   Mean   :  8223308   Mean   :66.69   Mean   :176.3  
##  3rd Qu.: 2.000   3rd Qu.: 10300000   3rd Qu.:70.00   3rd Qu.:200.0  
##  Max.   :21.000   Max.   :103000000   Max.   :77.00   Max.   :340.0
par(mfrow=c(2,4))
boxplot(bmi$erbmi, main="erbmi")
boxplot(bmi$ertpreat, main="ertpreat")
boxplot(bmi$ertseat, main="ertseat")
boxplot(bmi$euexfreq, main="euexfreq")
boxplot(bmi$eufastfdfrq, main="eufastfdfrq")
boxplot(bmi$eufinlwgt, main="eufinlwgt")
boxplot(bmi$euhgt, main="euhgt")
boxplot(bmi$euwgt, main="euwgt")

Vemos que todas las variables muestran valores atípicos. En el caso del BMI es esperable por su distribución. Otras variables, como ertseat tienen un número mucho mayor, por lo que debe tenerse en cuenta al modelizar y diagnosticar el modelo.

Una vez completada la descripción de las variables, convertimos a factor las variables categóricas:

factor.names <- c("erhhch", "erspemch", "ethgt", "etwgt", "eudietsoda", "eudrink", 
                  "eueat", "euexercise", "eufastfd", "euffyday", "eugroshp", 
                  "euinclvl", "eusnap",  "eumeat", "eumilk", "euprpmel", "eusoda", 
                  "eustores", "eustreason", "eutherm", "euwic",
                  "eeincome1", "erincome", "eufdsit", "eugenhth", "euincome2")
bmi[,factor.names] <- lapply(bmi[,factor.names], factor)

summary(bmi[,factor.names])
##  erhhch   erspemch  ethgt     etwgt     eudietsoda eudrink   eueat    
##  1: 500   -1:5249   0:10555   0:10535   -3:   1    -3:   1   -2:  55  
##  2: 213   1 : 221   1:   40   1:   53   -2:   4    -2:   7   1 :5822  
##  3:9924   2 :  90   2:   42   2:   49   -1:7767    1 :7164   2 :4760  
##           3 : 224                       1 :1121    2 :3465            
##           4 : 161                       2 :1677                       
##           5 :4692                       3 :  67                       
##                                                                       
##  euexercise eufastfd  euffyday  eugroshp  euinclvl eusnap    eumeat   
##  -2:   6    -2:  23   -2:   2   -3:   1   5:8760   -2:  35   -2:   9  
##  1 :6715    1 :6196   -1:4441   -2:   2   6:1877   1 :1091   -1:2973  
##  2 :3916    2 :4418   1 :2261   1 :6498            2 :9511   1 :6784  
##                       2 :3933   2 :2841                      2 : 871  
##                                 3 :1295                               
##                                                                       
##                                                                       
##  eumilk    euprpmel  eusoda    eustores    eustreason   eutherm   euwic    
##  -3:   1   -2:   3   -2:   3   -2:  50   2      :2895   -2:   3   -2:  22  
##  -2:   1   1 :6603   -1:3473   -1:2841   -1     :2891   -1:3853   -1:5130  
##  -1:2973   2 :2970   1 :2870   1 :5254   1      :2486   1 : 807   1 : 356  
##  1 : 150   3 :1061   2 :4291   2 :1927   3      :1030   2 :5974   2 :5129  
##  2 :7512                       3 : 337   4      : 668                      
##                                4 :  35   6      : 435                      
##                                5 : 193   (Other): 232                      
##  eeincome1 erincome  eufdsit   eugenhth  euincome2
##  -3: 102   -1: 213   -3:   1   -3:   2   -3: 197  
##  -2: 144   1 :6696   -2:  10   -2:  28   -2: 538  
##  1 :6696   2 : 504   1 :9987   1 :1943   -1:6577  
##  2 :3258   3 : 932   2 : 514   2 :3595   1 :1066  
##  3 : 437   4 :  33   3 : 125   3 :3288   2 :1917  
##            5 :2259             4 :1303   3 : 342  
##                                5 : 478

Búsqueda del modelo óptimo

Como hemos visto, el BMI se distribuye aproximadamente como una gamma, por lo que aplicamos este tipo de regresión. EM primer lugar, obtenemos modelos completos con todos los links disponibles, para luego hacer una selección de modelos.

full.id <- glm(erbmi ~ ., 
             data = bmi,
             family=Gamma(link="identity"))

par(mfrow=c(2,2)); plot(full.id); par(mfrow=c(1,1))

plot(residuals(full.id, type = "deviance"))

full.log <- glm(erbmi ~ ., 
             data = bmi, 
             family=Gamma(link="log"))

par(mfrow=c(2,2)); plot(full.log); par(mfrow=c(1,1))

plot(residuals(full.log, type = "deviance"))

full.inv <- glm(erbmi ~ ., # forumla
             data = bmi, # dataset
             family=Gamma(link="inverse"))

par(mfrow=c(2,2)); plot(full.inv); par(mfrow=c(1,1))

plot(residuals(full.inv, type = "deviance"))

Como vemos, los modelos completos cumplen condiciones como valores de los residuos DEVIANCE entre -2 y 2 en los tres link. Respecto al resto de diagnósticos, parece que el link identidad tiene mejores residuos.

Ahora seleccionamos modelos, para lo que hacemos una búsqueda mediante step. Usaremos el mejor modelo estimado por la función.

step.id.out <- step(full.id, direction = "backward")
step.log.out <- step(full.log, direction = "backward")
step.inv.out <- step(full.inv, direction = "backward")

Por tanto, los modelos que ajustaremos serán:

step.id.out$call
## glm(formula = erbmi ~ erincome + ethgt + etwgt + eusnap + eugenhth + 
##     eugroshp + euhgt + eusoda + eustreason + euwgt, family = Gamma(link = "identity"), 
##     data = bmi)
step.inv.out$call
## glm(formula = erbmi ~ eeincome1 + erhhch + erincome + erspemch + 
##     ethgt + etwgt + eueat + euexercise + euexfreq + eufastfd + 
##     eufinlwgt + eusnap + eugenhth + euhgt + euincome2 + eumeat + 
##     euprpmel + eusoda + euwgt + euwic, family = Gamma(link = "inverse"), 
##     data = bmi)
step.log.out$call
## glm(formula = erbmi ~ erspemch + ethgt + etwgt + eueat + euexfreq + 
##     eufinlwgt + eugenhth + euhgt + euincome2 + eusoda + euwgt + 
##     euwic, family = Gamma(link = "log"), data = bmi)
step.id <- glm(erbmi ~ erincome + ethgt + etwgt + eusnap + eugenhth + 
                 eugroshp + euhgt + eusoda + eustreason + euwgt, 
               family = Gamma(link = "identity"), 
               data = bmi)

par(mfrow=c(2,2)); plot(step.id); par(mfrow=c(1,1))

plot(residuals(step.id, type = "deviance"))

step.log <- glm(formula = erbmi ~ erspemch + ethgt + etwgt + eueat + euexfreq + 
                  eufinlwgt + eugenhth + euhgt + euincome2 + eusoda + euwgt + 
                  euwic, family = Gamma(link = "log"), data = bmi)
par(mfrow=c(2,2)); plot(step.log); par(mfrow=c(1,1))

plot(residuals(step.log, type = "deviance"))

step.inv <- glm(formula = erbmi ~ eeincome1 + erhhch + erincome + erspemch + 
                  ethgt + etwgt + eueat + euexercise + euexfreq + eufastfd + 
                  eufinlwgt + eusnap + eugenhth + euhgt + euincome2 + eumeat + 
                  euprpmel + eusoda + euwgt + euwic, family = Gamma(link = "inverse"), 
                data = bmi)
par(mfrow=c(2,2)); plot(step.inv); par(mfrow=c(1,1))

plot(residuals(step.inv, type = "deviance"))

En términos de AIC, el mejor modelo es el que emplea el enlace identidad, con AIC = 21647.03, contra el modelo con enlace logaritmo (AIC = 28302.22) y el de enlace inverso (AIC = 40643.39).

Si lo comparamos en cuanto a varianza explicada:

100 * (1 - (step.id$deviance/step.id$null.deviance))
## [1] 98.65191
100 * (1 - (step.log$deviance/step.log$null.deviance))
## [1] 97.47904
100 * (1 - (step.inv$deviance/step.inv$null.deviance))
## [1] 91.98856

De nuevo, el mejor modelo es el que emplea en enlace identidad, que explica un 98.65% de la DEVIANCE. Si observamos sus coeficientes:

summary(step.id)
## 
## Call:
## glm(formula = erbmi ~ erincome + ethgt + etwgt + eusnap + eugenhth + 
##     eugroshp + euhgt + eusoda + eustreason + euwgt, family = Gamma(link = "identity"), 
##     data = bmi)
## 
## Deviance Residuals: 
##       Min         1Q     Median         3Q        Max  
## -0.131176  -0.010936   0.001528   0.009796   0.229545  
## 
## Coefficients:
##                Estimate Std. Error  t value Pr(>|t|)    
## (Intercept)  51.4083315  0.8729159   58.893  < 2e-16 ***
## erincome1     0.0745781  0.0443039    1.683 0.092341 .  
## erincome2     0.1312195  0.0521158    2.518 0.011822 *  
## erincome3     0.0839181  0.0484358    1.733 0.083203 .  
## erincome4     0.1090085  0.1143751    0.953 0.340572    
## erincome5     0.1567179  0.0461082    3.399 0.000679 ***
## ethgt1        0.8722284  0.0924147    9.438  < 2e-16 ***
## ethgt2        0.9540602  0.1216994    7.839 4.96e-15 ***
## etwgt1       -1.4951243  0.1708827   -8.749  < 2e-16 ***
## etwgt2       -0.9310082  0.0646855  -14.393  < 2e-16 ***
## eusnap1       0.0017920  0.1044026    0.017 0.986306    
## eusnap2      -0.0698268  0.1021802   -0.683 0.494389    
## eugenhth-2    0.2271269  0.4625373    0.491 0.623404    
## eugenhth1     0.3133602  0.4455386    0.703 0.481866    
## eugenhth2     0.3217543  0.4454499    0.722 0.470118    
## eugenhth3     0.3542196  0.4454281    0.795 0.426495    
## eugenhth4     0.3976043  0.4456215    0.892 0.372281    
## eugenhth5     0.4523042  0.4461883    1.014 0.310746    
## eugroshp-2    0.5676538  1.0056308    0.564 0.572443    
## eugroshp1     0.7083560  0.8848995    0.800 0.423443    
## eugroshp2     0.7115029  0.8897742    0.800 0.423935    
## eugroshp3     0.6458298  0.8848067    0.730 0.465461    
## euhgt        -0.7684382  0.0019201 -400.199  < 2e-16 ***
## eusoda-1     -0.1218794  0.3278954   -0.372 0.710122    
## eusoda1      -0.1318953  0.3279686   -0.402 0.687576    
## eusoda2      -0.1569894  0.3278840   -0.479 0.632094    
## eustreason-2 -0.7734256  0.5980125   -1.293 0.195926    
## eustreason-1 -0.8762340  0.5992355   -1.462 0.143701    
## eustreason1  -0.8137849  0.5921778   -1.374 0.169401    
## eustreason2  -0.8602364  0.5921865   -1.453 0.146352    
## eustreason3  -0.8598634  0.5924067   -1.451 0.146677    
## eustreason4  -0.7647428  0.5926035   -1.290 0.196912    
## eustreason5  -0.7910744  0.5942001   -1.331 0.183110    
## eustreason6  -0.8720373  0.5928597   -1.471 0.141348    
## euwgt         0.1560453  0.0002035  766.798  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Gamma family taken to be 0.0006107502)
## 
##     Null deviance: 476.2832  on 10636  degrees of freedom
## Residual deviance:   6.4207  on 10602  degrees of freedom
## AIC: 21647
## 
## Number of Fisher Scoring iterations: 5

Vemos que solo los coeficientes correspondientes a las variables erincome, ethgt, euhgt y euwgt, que se refieren a ingresos, peso y altura, son significativos.

Al emplear un enlace identidad, podemos interpretar el modelo como que cada incremento de una unidad en \(X_j\) hace crecer el valor esperado del BMI, \(E(BMI)\) en \(\beta_j\), asumiendo que el resto de variables se mantienen constantes.

Comparación con un modelo lineal

A continuación, ajustaremos el mismo modelo que hemos estimado como óptimo para ver si hay diferencias.

gaussian <- glm(erbmi ~ erincome + ethgt + etwgt + eusnap + eugenhth + 
                 eugroshp + euhgt + eusoda + eustreason + euwgt,
             family = ("gaussian"),
             data = bmi)

par(mfrow=c(2,2)); plot(gaussian); par(mfrow=c(1,1))

Los residuos no se ajustan a la normalidad, además de ser heterocedásticos. Por tanto, vemos que la elección de un modelo gamma es más adecuado a la hora de describir los datos.

Validación cruzada

Pese a que el objetivo del modelo no es la predicción, realizamos una verificación mediante validación cruzada para comprobar como se ajustan los valores.

set.seed(1)
data_partition <- createDataPartition(bmi$erbmi, 
                                      times = 1, p = 0.8, 
                                      list = FALSE)

train <- bmi[data_partition,]
test  <- bmi[-data_partition,]

cv.model <- glm(erbmi ~ erincome + ethgt + etwgt + eusnap + eugenhth + 
                  eugroshp + euhgt + eusoda + eustreason + euwgt, 
                family = Gamma(link = "identity"), 
                data = test)

prediction <- predict(cv.model, newdata=test, type="response")
plot(test$erbmi, prediction)

Vemos que la relación es prácticamente lineal, aunque los valores extremos tienen más variabilidad.